In [ ]:
import matplotlib.pyplot as plt
import cobra
from cobra.io import validate_sbml_model
import importlib
from xml.etree import ElementTree
import utils.Model_correction as mc
import sys
import utils.model_maj as mj
import utils.viz_utils as vu
#import cplex
import cbmpy
import pandas as pd
import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
from tqdm import tqdm
from itertools import cycle
import seaborn as sns

INFO: No xlwt module available, Excel spreadsheet creation disabled


CBGLPK based on swiglpk: not all methods implimented yet!

*****
Using CPLEX
*****


INFO: No xlrd module available, Excel spreadsheet reading disabled


***********************************************************************
* Welcome to CBMPy (0.8.2) - PySCeS Constraint Based Modelling        *
*                http://cbmpy.sourceforge.net                         *
* Copyright(C) Brett G. Olivier 2014 - 2020                           *
* Systems Biology Lab, Vrije Universiteit Amsterdam                   *
* Amsterdam, The Netherlands                                          *
* CBMPy is developed as part of the BeBasic MetaToolKit Project       *
* Distributed under the GNU GPL v 3.0 licence, see                    *
* LICENCE (supplied with this release) for details                    *
***********************************************************************
In [ ]:
importlib.reload(vu)
importlib.reload(mc)
Out[ ]:
<module 'utils.Model_correction' from '/home/jerem/Documents/master/RBA_cutsets/code/utils/Model_correction.py'>

Models loading and bounds modification¶

In [ ]:
HepG2, errors = validate_sbml_model("../models_storage/HepG2_v8.xml")
In [ ]:
iHep, errors = validate_sbml_model("../models_storage/iHep_v3_genenames.xml")
In [ ]:
HepG2_C = HepG2.copy()
In [ ]:
HepG2.reactions.EX_m01965x.bounds = (-1000.0,-1000.0) #Glucose exchange, set to the maximal input value observed with FVA
#iHep.reactions.EX_m01965x.bounds = (-1000.0, -1000.0)
iHep.reactions.HMR_4281.bounds = (0.0,0.0) # Blocking peroxysome fermentation in healthy model
In [ ]:
HepG2.reactions.HMR_3883.bounds = (0.0,0.0) #Serine --> Pyruvate exchange
HepG2.reactions.EX_m02630x.bounds = (0.0,0.0) #O2 exchange
HepG2.reactions.EX_m02819x.bounds = (0.0,1000.0) #Pyruvate exchange
HepG2.reactions.EX_m01910x.bounds = (0.0,0.0) #Galactose exchange
HepG2.reactions.EX_m01743x.bounds = (0.0, 1000.0) #D-Ribulose exchange
HepG2.reactions.EX_m01962x.bounds = (0.0,1000.0) #glucosamine exchange
HepG2.reactions.HMR_4316.bounds = (-1000.0,1000.0) # Glucose --> D-Glucitol
HepG2.reactions.EX_m02896x.bounds = (0.0,1000.0) # Serine intake
HepG2.reactions.EX_m01682x.bounds = (-1000.0,0.0) # Glucitol secretion blocked. Kept the intake just in case.
HepG2.reactions.EX_m01840x.bounds = (-1000.0,0.0) #Fructose exchange
HepG2.reactions.HMR_4300.bounds = (494.6155049539209,494.6155049539209)
#HepG2.reactions.HMR_4281.bounds = (100.0,1000.0) #Fermentation perox.
#HepG2.reactions.HMR_4930.bounds = (-1000.0,1000.0) # Pyruvate transfer from cytoplasm to peroxysome
#HepG2.reactions.HMR_4281.bounds = (0.0,0.0) # Fermentation in peroxysome

Optimisation¶

In [ ]:
iHep.solver = "cplex"
iHep.objective = "biomass_components"
sol_iHep = iHep.optimize()
HepG2.solver = "cplex"
HepG2.objective = "biomass_components"
sol_G2 = HepG2.optimize()
HepG2_C.solver = "cplex"
HepG2_C.objective = "biomass_components"
sol_G2_C= HepG2_C.optimize()

Heatmap¶

In [ ]:
df_Ci, df_Ri, df_Si, df_Mi, df_Pi, df_Xi, df_Li, df_Gi, df_Ni = vu.build_reaction_df(iHep)
df_C, df_R, df_S, df_M, df_P, df_X, df_L, df_G, df_N = vu.build_reaction_df(HepG2)
subS_HepG2 = vu.get_subsystem_fluxes(vu.build_reaction_df(HepG2))
subS_iHep = vu.get_subsystem_fluxes(vu.build_reaction_df(iHep))
subS_G2C = vu.get_subsystem_fluxes(vu.build_reaction_df(HepG2_C))
df_both = pd.concat([subS_HepG2, subS_G2C, subS_iHep], axis = 1)
df_both.columns = ["HepG2","HepG2_copy", "iHep"]
#df_both = df_both.loc[((df_both["iHep"]/df_both["HepG2"] >= 2.0) | (df_both["HepG2"]/df_both["iHep"] >= 2.0)) & ((df_both["HepG2"]>=1000.0) | (df_both["iHep"]>=1000.0))]  # --> Filtre relatif + absolu 
df_both = df_both.loc[(df_both["HepG2"]>=1000.0) | (df_both["iHep"]>=1000.0)] # --> Filtre absolu
clustermap = sns.clustermap(df_both.fillna(0), figsize =(10,20))

Heatmap : Représentation sous forme d'une heatmap des flux de chaque sous-systèmes actif dans trois modèles de cellule du foie différentes :

  • iHep : Modèle sain
  • HepG2 : modèle cancéreux, aux réactions corrigées
  • HepG2_copy : modèle cancéreux, aux réactions non corrigées.

On s'attend à observer une distribution des flux dirigée vers la production d'énergie, de nucléotides, d'acides aminés ou d'acides gras pour les modèles cancéreux.

La glycolyse est effectivement plus active dans le modèle HepG2 corrigé. \ On observe également une activité importante du métabolisme du folate dans le modèle HepG2 (métabolisme des nucléotides)\ Cependant, l'activité du métabolisme des pentoses phosphates reste bien trop faible, il en va de même pour le métabolisme du pyruvate.

Comparaison, par compartiment, des subsystems les plus actifs pour les deux modèles :¶

On s'attend à observer les différences suivantes dans le modèle cancéreux par rapport au modèle sain :

  • activité accrue du métabolisme des nucléotides
  • activité accrue de la voie des pentoses phosphates
  • activité accrue de la glycolyse
  • activité accrue de la voie de biosynthèse des hexosamines
  • activité accrue de la lipogénèse

--> Résultats :

  • Activité deux fois plus importante des réactions de glycolyse dans le cytoplasme du modèle cancéreux, par rapport au modèle sain.
  • Activité bien plus importante de la voie des pentoses phosphates et du métabolisme des nucléotides dans le modèle sain.
  • Activité 2x plus importante de l'activité du métabolisme de l'acide folique dans le cytoplasme du modèle cancéreux.
  • Aucune différence significative repérée pour la lipogénèse ou la biosynthèse des hexosamines

Comparaison des flux attribués à chaque groupe de réactions entre le modèle sain et le modèle cancéreux

In [ ]:
r = HepG2.reactions.HMR_0162
r.subsystem
Out[ ]:
'Carnitine shuttle (mitochondrial)'

Barplots comparaison par groupe des flux associés à chaque gène¶

In [ ]:
barplots_dict = vu.subsystem_barplots(iHep, HepG2_C, 'iHep', 'HepG2_copy',{'glycolysis / gluconeogenesis' : (10,14), 'pyruvate metabolism' : (7,11), 'pentose phosphate pathway' : (5,7), 'folate metabolism' : (5,7), 'carnitine shuttle' : (5,7), 'beta oxidation' : (5,7), "nucleotide metabolism" : (10,14)})
In [ ]:
barplots_dict_2 = vu.subsystem_barplots(iHep, HepG2, 'iHep', 'HepG2',{'glycolysis / gluconeogenesis' : (10,14), 'pyruvate metabolism' : (7,11), 'pentose phosphate pathway' : (5,7), 'folate metabolism' : (5,7), 'carnitine shuttle' : (5,7), 'beta oxidation' : (5,7), "nucleotide metabolism" : (10,14)})

Nucleotide metabolism¶

On s'attend globalement à une activité plus importante du métabolisme des nucléotides dans le modèle cancéreux. Cependant, on n'observe une activité plus importante des protéines cancéreuses que lorsqu'elles sont impliquées dans la synthèse d'ADP.

Cela pourrait indiquer un potentiel biais : l'ADP est indiqué comme produit dans la fonction de biomasse données en objectif au modèle cancéreux.

In [ ]:
barplots_dict_2["nucleotide metabolism"][0]
Out[ ]:
In [ ]:
barplots_dict["nucleotide metabolism"][0]

Glycolysis reactions¶

--> On s'attend à observer :

  • une activité très réduite de FBP1, PCK1/2, PKLR, GPI, PFKL, ALDOA/B/C, TPI, PC
  • une activité similaire entre modèle cancéreux et sain de GAPDH, PGK1/2, PGAM1, ENO1/2/3, LDHA.

--> On observe :

  • Correct : Activité importante de la pyruvate kinase chez les deux types cellulaires : gènes PKLR/PKM.
    • en théorie, PKLR est prédominante dans les cellules saines, et PKM dans les cellules cancéreuses. La nature booléenne de l'association GPR ne permet pas de différencier les deux.
  • Correct : activité légèrement plus importante de TPI Chez le modèle sain par rapport au modèle cancéreux.
  • Correct : Activité globalement plus importante pour les enzymes de la glycolyse dans le modèle cancéreux :
    • GPI, HK1/2/3 --> + de glycolyse et donc plus de fermentation, et d'activité vers la voie des pentoses phosphates.
In [ ]:
barplots_dict_2["glycolysis / gluconeogenesis"][0]
Out[ ]:
In [ ]:
barplots_dict['glycolysis / gluconeogenesis'][0]

Pyruvate metabolism :¶

Cette figure sert surtout en conjonction avec la précédente où on peut voir que les flux associés aux gènes PKLM et PKR sont les mêmes entre le modèle sain et cancéreux.

In [ ]:
barplots_dict_2["pyruvate metabolism"][0]
Out[ ]:
In [ ]:
barplots_dict["pyruvate metabolism"][0]
Out[ ]:

Pentose phosphate pathway¶

On s'attend à une activité de la voie des pentoses phosphates plus importante chez la cellule cancéreuse que chez la cellule saine, pourtant ce n'est pas du tout le cas.

La seule activité plus importante pour le modèle cancéreuse est celle de la phosphofructokinase. Or il s'agit d'une enzyme censée être sous-exprimée dans une cellule cancéreuse.

In [ ]:
barplots_dict_2["pentose phosphate pathway"][0]
Out[ ]:
In [ ]:
barplots_dict["pentose phosphate pathway"][0]
Out[ ]:

Folate metabolism¶

Le métabolisme du folate est supposé être suractivé dans une cellule cancéreuse.

C'est bien ce qu'on observe ici avec une activité beaucoup plus importante de DHFR, qui est d'ailleurs une cible thérapeutique anticancéreuse connue.

In [ ]:
barplots_dict_2["folate metabolism"][0]
Out[ ]:
In [ ]:
barplots_dict["folate metabolism"][0]
Out[ ]:

Beta oxidation¶

On s'attend à une activité très réduite des protéines de la beta-oxydation (ici, toutes les beta oxydations sont regroupées : mitochondriale et peroxysomale)

On observe effectivement une activité inexistante (flux nuls)

In [ ]:
barplots_dict_2["beta oxidation"][0]
Out[ ]:
In [ ]:
barplots_dict["beta oxidation"][0]
Out[ ]:

Carnitine shuttle¶

On s'attend à une activité très réduite des protéines associées à la navette carnitine.

On observe effectivement une activité très réduite de ces protéines.

In [ ]:
barplots_dict_2["carnitine shuttle"][0]
Out[ ]:
In [ ]:
barplots_dict["carnitine shuttle"][0]
Out[ ]:

Treemaps -- BY SUBSYSTEM¶

In [ ]:
Subsystem_fluxes_cancer = vu.build_reaction_df(HepG2, False)
treemap_dataframe = pd.DataFrame(Subsystem_fluxes_cancer["Glycolysis / Gluconeogenesis"])
title = "Glycolysis reactions -- <b>CANCER MODEL<br />"

vu.plot_treemap(treemap_dataframe, HepG2, title, path = ["compartment", "id"], color_by="compartment")
In [ ]:
Subsystem_fluxes_healthy = vu.build_reaction_df(iHep, False)
treemap_dataframe = pd.DataFrame(Subsystem_fluxes_healthy["Glycolysis / Gluconeogenesis"])
title = "Glycolysis reactions -- <b>HEALTHY MODEL<br />"

vu.plot_treemap(treemap_dataframe, iHep, title, path = ["compartment", "id"], color_by="compartment")

Treemaps + barplots -- BY COMPARTMENT¶

Cytoplasme¶

In [ ]:
title = "Cytoplasm reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_C,HepG2,title).show(renderer = 'notebook')
In [ ]:
title = "Cytoplasm reactions -- <b>HEALTHY MODEL <br />"
vu.plot_treemap(df_Ci, iHep, title).show(renderer = "notebook")
In [ ]:
barplots_dict["C_c"][0]
Out[ ]:

Mitochondrie¶

In [ ]:
title="Mitochondrial reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(df_Mi, iHep, title).show(renderer = "notebook")
In [ ]:
title="Mitochondrial reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_M, HepG2, title).show(renderer = "notebook")
In [ ]:
barplots_dict["C_m"][0]
Out[ ]:

Peroxysome¶

In [ ]:
title="Peroxysomal reactions -- <b>HEALTHY MODEL<br />"
vu.plot_treemap(df_Pi, iHep, title).show(renderer = "notebook")
In [ ]:
title="Peroxysomal reactions -- <b>CANCER MODEL<br />"
vu.plot_treemap(df_P, HepG2, title).show(renderer = "notebook")
In [ ]:
barplots_dict["C_p"][0]
Out[ ]: